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Reconstruction of the dynamics (quantum process tomography) of the single-exciton manifold in 
energy transfer systems is proposed here on the basis of two-dimensional fluorescence spectroscopy 
(2D-FS) with phase-modulation. The quantum-process-tomography protocol introduced here bene- 
hts from, e.g., the sensitivity enhancement ascribed to 2D-FS. Although the isotropically averaged 
spectroscopic signals depend on the quantum yield parameter F of the doubly-excited-exciton mani¬ 
fold, it is shown that the reconstruction of the dynamics is insensitive to this parameter. Applications 
to foundational and applied problems, as well as further extensions, are discussed. 
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I. INTRODUCTION 

Since Chuang and Nielsen’s 1996 seminal proposal to 
experimentally reconstruct the evolution operator of a 
quantum black box [1], a variety of experiments have 
been proposed [2-4] and some of these have been imple¬ 
mented [5-7] . Although most of these ideas involved rela¬ 
tively “clean” optical systems in the context of quantum 
information processing (QIP), recently there has been 
growing interest in the application of quantum process to¬ 
mography to study electronically coupled molecular sys¬ 
tems [8-10]. Such combined theoretical and experimental 
studies on excitation energy transfer serves to bring to¬ 
gether the QIP and physical chemistry communities. 

The current version of QPT for excitonic systems 
[8, 9, 11] is based on the method of 2D Photon Echo 
Spectroscopy (2D-PES) [12, 13]. It therefore relies on 
the wave-vector phase-matching condition, which works 
for macroscopic systems with many chromophores (see 
Refs. [11-13] for details). For this reason, the proposal de¬ 
veloped in Refs. [8-10] is not suitable for single-molecule 
QPT, and it is hence desirable to implement the phase¬ 
cycling (PCT) or phase-modulation techniques (PMTs) 
[14, 15]. In most 2D-PES experiments, the system inter¬ 
acts with three non-collinear ultrafast laser pulses, which 
gives rise to a third-order polarization that propagates 
in the wave-vector matched direction [11-13]. The trans¬ 
mitted signal must be separated from background laser 
light, which is inadvertently scattered by the sample in 
the same direction as the third-order signal. The pres¬ 
ence of background scattering is a limiting factor to the 
sensitivity of 2D-PES experiments. In 2D-FS, the system 
interacts with four collinear laser pulses, and the signal is 
detected by monitoring nonlinear contributions to the en¬ 
suing fluorescence signal [15, 16]. The red-shifted fluores¬ 
cence can be easily separated from background scattered 
laser light by using long-pass spectral filters. Specific non¬ 
linear contributions to the fluorescence signal are isolated 
according to the phase-modulation schemes described in 


previous work [15, 16]. Because the 2D-FS method is 
based on the detection of incoherent fluorescence signals, 
it may be applied to systems of small numbers of chro¬ 
mophores, quantum dots, or thin film materials. 

Despite the high-sensitivity advantages afforded to 
fluorescence-detected PMTs, which are useful for stud¬ 
ies of biological molecules and molecular aggregates [17- 
19], these methods have been less commonly practiced 
than four-wave mixing approaches to 2D-PES. However, 
recent theoretical [15] and experimental [16, 20] progress 
with classical light have enabled PMTs for a variety of 
complex molecular systems relevant to exciton dynam¬ 
ics. These recent developments, and the general theory 
of open quantum systems- i.e., quantum systems coupled 
to the environment- are combined here to formulate a 
self-consistent theory of QPT that is based on collinear 
PMT with synchronous detection. 

For the ideal situation when nonradiative processes 
are neglected in the doubly-excited-exciton manifold, the 
quantum yield parameter of this manifold is set to P = 2. 
Under this circumstance, it was shown that 2D-FS coin¬ 
cides with 2D-PES [15, 16, 20]. It is shown below that 
this equivalence also holds at the level of quantum pro¬ 
cess tomography, i.e., the protocol introduced here gen¬ 
eralizes the protocol in Refs. [8, 9] to the more realistic 
situation 0 < P <2. 


II. INITIAL CONSIDERATIONS ON QPT, 
SYSTEM MODEL AND 2D-FS 

Before introducing the reconstruction of the dynam¬ 
ics, it is necessary to state some remarks on the basics 
of process tomography, the system-of-interest model and 
2D-FS. 

Quantum Process Tomography Tensor —In quantum 
mechanics, the state of a physicochemical system S is 
described by a density operator p. Time evolution of 
quantum states is governed by the Schrodinger equation. 
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which is linear in the state of the system. This linearity 
allows for a description of the system’s dynamics in terms 
of a linear map, Xt ■ Po Pt- After projecting onto a 
complete orthonormal basis {|n)}, the map reads 

{n\p{t)\m) = Xnmi.^^{t){v\p{0)\p), ( 1 ) 

where x-nmi^nit) stands for the process tomography 
tensor [2, 8, 21]. For Hamiltonian dynamics with 
H\n) = A„|n), Xnm,.^.{t) = [22- 

24]. Thus, population-to-coherence [xurni^i^it)] and the 
reverse [xnni^fj,it)] process are prevented by the Kronecker 
deltas Sni/Smfj,- Clearly, this restriction is not present if 
driving fields are present or if the system of interest is 
coupled to its environment [22-24]. 

In the general case of open quantum systems, the func¬ 
tional form of Eq. (1) remains valid under some condi¬ 
tions. (i) If the coupling to the bath is weak, Eq. (1) 
holds for Markovian and non-Markovian processes and 
the process tensor is independent of the initial state 
(see, e.g., Refs. [8, 9] and references therein), (ii) If the 
coupling to the bath is strong, and the initial system- 
environment correlations cannot be neglected, Eq. (1) 
holds after including those initial correlations in Xnmuft {t) 
(see Refs. [22-24] for details), (iii) Because the initial 
bath correlations vanish at high temperature, even for 
strong coupling [25], then Xnmvp.{t) can be defined in¬ 
dependently of the initial state in the strong coupling 
regime entered at high temperatures. 

After identifying the conditions under which Eq. (I) 
holds, it is relevant to consider some of the main proper¬ 
ties of the QPT tensor [8], namely. 


* 

Xmniip, Xmnp.i'i 

(2) 

'y ' Xnnp.i'iT') = Spu, 
n 

(3) 


(4) 


nmi'^ 


where z is any complex valued vector. Equation (2) en¬ 
sures the Hermitian character of the density operator, 
p = p\ while Eq. (3) guaranties probability conservation, 
tr/3(t) = I. The last property is a consequence of the fact 
that p{t) remains positive-semidefinite under unitary op¬ 
erations. 

The objective of QPT is the experimental reconstruc¬ 
tion of the process tomography tensor Xnmvii{t). 

Model —Consider an excitonic dimer described by Hg 
and given by 

Hq = widjdi -b u}2ala2 + J {^ 1^2 + d\a^ , (5) 

where aj and hi are the creation and annihilation op¬ 
erators for site i, vji yf VJ 2 are the site energies while 
J ^ 0 is the Coulombic coupling between chromophores. 
By defining the average frequency w = ^{wi + ^ 2 ), 


the half-difference A = — W 2 ) and the mixing 

angle B = ^ arctan( J/A), it is possible to introduce 
the creation and annihilation operators, Cp = cos0ai -I- 
sin 0 a 2 and = sin^dj + cos 9a\, of the p-th delocal¬ 
ized exciton state with energy Wp = ro ± A sec 29 and 
p G {e,e'}. Starting from the ground state \g), the 
single-exciton states are conveniently defined as | e) = 
cl\g) and | e') = Cg, | 5 ), while the biexciton state as 
]/) = a\al\g) = 4cj],|p) with zuf = -b 072 = 
uje + vUe'- The dipole vectors at each site are set to 
di = diBz and d 2 = d2COs{4>)ez + c ?2 sin((/))ea;. So 
that, pL^g = c ?2 sin 9 sin (pe^ + (c?i cos 0 -b c ?2 sin 9 cos (f) , 

fj-e'g = ^2 COS^sinc^Ga; -b ( —di sin0-b ^2 COS0COS(^) e^, 

M/e = ^2 COS 9 sin ipBx + {di sin 9 + d 2 cos 9 cos (j)) and 
fife' = —d 2 sin 9 sin -b (di cos 9 — d 2 sin 9 cos (f)) . 

Although exciton-exciton binding or repulsion terms 
are not included here, it is considered that each excitonic 
manifold contributes to the spectroscopic signal with a 
weight given by their fluorescence quantum yield coef¬ 
ficients r,y. Specifically, it is assumed that the quan¬ 
tum yield of the two singly excitonic states are the same 
and equal to I, while for the doubly excitonic manifold, 
it is assumed that P/ = P with 0 < P < 2. In the 
ideal case in which two photons are emitted via the path 
I /) —>• I e, e') —>■ I p), P = 2. It is possible that singlet- 
singlet annihilation would convert a doubly-excited ex- 
cition into a singly-excited exciton [26], which in the 
absence of non-radiative decay would result in P = I. 
However, because of the abundance of non-radiative re¬ 
laxation pathways for highly excited states, the quantum 
yield of the doubly-excitonic manifold is expected to be 
smaller than that of the singly excitonic manifold, so that 
values smaller than unity are expected. For example, for 
membrane-supported self-assembled porphyrin dimers, it 
was found that P = 0.31 [16, 20]. 

For convenience, the dimer Hamiltonian can be writ¬ 
ten as Hs = J2,y={g e e' \^){ ^1- To account for the 
influence of the local vibrational environment in the exci¬ 
tonic dimer, coupling to a thermally equilibrated phonon 
bath at inverse temperature j3 is considered next. Specif¬ 
ically, the Hamiltonian of the environment is given by 

He = Ep=e,e' J2n ^n,p + l/s), where de¬ 

notes the frequency of the environment modes. The in¬ 
teraction is described by^sE = Ae|e)(e|-b Ee' |e')(e'|-b 
(^Ee + Ee'^ \e'){e'\ with Ep = (^«.p + ■ 

p and bn,p are the creation and annihilation bosonic op¬ 
erators of the n—th mode of the vibrational environment 
in the p—site. Xn,p measures the interaction strength be¬ 
tween the n—th mode of the environment and the p—th 
site. The net effect of the local environment is encoded 
in the spectral density J„ = ^n,p^n,p^(.^ ~ ^n)- 

2D Fluorescence Spectroscopy (2D-FS) —The main dif¬ 
ference between the QPT scheme introduced below and 
previous QPT proposals is the spectroscopic technique, 
2D-FS, which the present proposal is based on. It is 
therefore relevant to discuss the main differences and ad- 
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vantages that 2D-FS has over, e.g., 2D-PES. The 2D-FS 
method isolates the nonlinear optical response of a mate¬ 
rial system by monitoring fluorescence signals. Because 
fluorescence can be efficiently separated from background 
scattered light, the 2D-FS approach can be used to per¬ 
form experiments that require very high detection sensi¬ 
tivity [15, 16, 20]. Moreover, the collinear beam geometry 
used in 2D-FS has the advantage that every illuminated 
molecule experiences the same optical phase condition 
at every instant in time. Since the incoherent fluores¬ 
cence signal is emitted isotropically, very small numbers 
of molecules may be studied in this way. [15, 16, 20]. 

The 2D-FS observable is proportional to the fourth- 
order excited populations, 

(i(t)) = trip(4)(t), (6) 

with A = X]i/={e e' /} generated by the action 

of the operator V(t') that comprises the excitation by 
four weak non-overlapping laser pulses. 


V(t') 


4 

—A p,-eiE(t' — t) 

i=l 




• (7) 


Here A denotes the maximum intensity of the pulses’ elec¬ 
tric field, and p, the dipole operator, e^, ti, uJi and (pi 
stand for the polarization vector, time center, frequency 
and phase of the i-th laser pulse. The pulse envelope 
E{t) is chosen to be Gaussian with fixed width tr, i.e., 
E{t) = e“‘ . In the model under consideration, the 

only optically allowed transitions are between states dif¬ 
fering by one excitation. Hence, the only non-vanishing 
dipole transition matrix elements are jiij = fiji with 
= {^ 5 , e'g,/e,/e'}. Details about the derivation and 
the explicit functional form the fourth-order density ma¬ 
trix can be found in Appendices A and B. 

For the purpose of extracting the QPT tensor from 
the 2D-FS experimental signals, only the rephasing sig¬ 
nals with global phase t^reph = —4>i + 4>2 + < 1)3 ~ 4^4 
be considered below (see Fig. 1). Thus, assuming that 
the rotating wave approximation (RWA) holds, the inter¬ 
actions with the electromagnetic fields are characterized 
by 

Vi=-XP< ( 8 ) 

V2 = -Xp> ■ e2E{t - (9) 

Fa = -AA> ■ esEit - ( 10 ) 

V 4 = -AA< ■ e4E(t - (11) 

where = J2Ljp<Ljq t^pq\p){Q\ promotes emissions from 
the ket and absorptions on the bra, and p^ = {p'^Y in¬ 
duces the opposite processes. For this particular selection 
of the global phase, t^reph = + ^2 + 4‘3~ the 2D- 

FS signals are equivalent to the rephasing spectroscopic 
signals in the photon-echo direction kpg = —ki -|-k 2 Aka 
when F = 2 [15, 16, 20]. It is shown below that 


the present QPT protocol reduces to the protocol in 
Refs. [8, 9] when F = 2 as well. 


III. 2D-FS QPT 

As stated above, the main goal of QPT is the recon¬ 
struction of the dynamics of the density operator. In do¬ 
ing so, it is assumed that the structural parameters of the 
model, namely, the transition frequencies Wij = Wi — Wj 
and the electric dipole transition matrix element jjiij are 
all known. This prerequisite is not an issue because in¬ 
formation about the transition frequencies is routinely 
obtained from linear absorption spectra, and the tran¬ 
sition dipole directions can be inferred form structural 
measurements and polarization spectroscopy [ 8 ] . 

Once the structural parameters are defined, the recon¬ 
struction of the dynamics comprises three main parts: 
(i) initial state preparation, (ii) evolution and (iii) final 
state detection. In describing these stages, it is useful to 
introduce the standard time intervals {r, T, t} instead of 
the time center ti of each pulse [12, 13]. The time differ¬ 
ence between the second and the first pulse defines the 
coherence time interval r = t 2 — ti. The time interval 
between the third and the second pulse, T = t^ — t 2 is 
known as the waiting time, which defines the quantum 
channel to be characterized by the QPT scheme. Finally, 
the difference between the fourth and the third pulse, 
t = ti — ts, denotes the echo time. 

Initial State Preparation —The excitonic system, be¬ 
fore any electromagnetic perturbation, is assumed to be 
in the ground state, /5(— 00 ) = \g){g\. Thus, the basic 
idea is to make use of the first two pulses to prepare the 
effective initial density matrix at T = 0, = 0), 

and use the last two pulses to read out the state. 

After applying second order perturbation theory in A, 
and under the assumption that the RWA holds in this 
case (see Appendix A for details), the effective initial 
state reads 

E C'£iC'«,(Mpg-ei)(M,g-e2) 

p,gG{e,e'} ( 12 ) 

X Ggpi.T){\q){p\ -5pq\g){g\), 

where Qijij) is the propagator of the optical coherence 
\i){j\. For simplicity, it can be assumed as Qijir) = 
0(t) exp[(—iwij — Fij)T] begin Fy dephasing rates, and 
the Heaviside function 0(r) ensures causality. The co¬ 
efficients are purely imaginary and given by CJ. = 

iXV2TT^e~'^ Because at this level there is no 

influence of the doubly-excited exciton manifold, the ef¬ 
fective initial state in Eq. (12) coincides with the effective 
initial state prepared by 2D-PES in Ref. [ 8 ]. 


The contributions to in Eq. (12) are clearly 

depicted in the lower panel of Fig. 1 (see Figs. A, B, C 
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Final state detection 





Initial state preparation 

/^e r^e' 

I. 1-1 ^/.In 


\ 9){a \^4>2 1 e )( e 1 A 

B le'Xel 

l9)(e| 02^5 )(e| 

l.9>(fil'^0l l.9)(!?i^0i 

02 1 5)( e 1 

9 X a 1^01 


r^e r^e 


C^e. /^e 

^CJi ^(jJ2 
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D leXe'l 
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FIG. 1. Double-sided Feynman’s diagrams for the initial state preparation (lower panel) and state detection (upper panel) 
that lead to the rephasing spectroscopic signals [S'Fs]e(,’^ 2 ^,’il^% 4 ’‘(r, T, t) in Eq. (14). The 2D-FS signals are synchronously-phase 
detected with respect to the modulated laser fields at frequency drep = —di + <(>2 + da — 4’i- la the lower and npper panels, the 
diagrams are grouped according to the probability they occur and respectively. 


and D). Starting from the system the ground state |g)(g|, 
in the RWA and selecting only those contributions with 
—-b (/> 2 , the first pulse can only excite the bra and then 
creates an the optical coherence |5)(p| with probability 
CP^. This coherence \g){p\ evolves under the action of 
Qgp{T) for a time r when the second pulse prepares the 
state |g)(p| with probability or a hole — \g){g\ with 
probability . 

As in the case of QPT based on 2D-PES [8, 9], to pre¬ 
pare the set of four linearly independent states in Eq. (12) 
(see also Figs. l.A-l.D), it suffices to consider a pulse tool¬ 
box of two waveforms with carrier frequencies {a;+,aj_} 
that create | e) and | e') with different amplitudes. Of 
course, the discrimination in the preparation of | e) and 
I e') depends on how close to resonance the carrier fre¬ 
quencies are. For an extensive and detailed analysis on 
this respect, see Refs. [8, 9]. 

Evolution —Once the initial state is effec¬ 

tively prepared, i.e., after the action of the first-two 
pulses T > 3(7, the system evolves over a time T un¬ 
der the action of the super operator xiT), according to 

/3e^’e|(r) = x(T)p^^’^2^(0). (13) 

To avoid contamination of the initial state by terms pro¬ 
portional to a hole every time there is a single-exciton 
population \p){p\, it is assumed that {ab\x(T)\gg) = 


Xabgg{T) = SagSbg, which is equivalent to neglect pro¬ 
cesses where phonons induce upward optical transitions 
and spontaneous excitation from the single to the dou¬ 
ble exciton manifolds [8]. Up to this condition, x{T) in 
Eq. (13) describes the dynamics induced by any bath 
model and accounts for any system-bath coupling. 

Final State Detection —The very nature of the fluores¬ 
cence detection in 2D-FS suggests considering contribu¬ 
tions from excitation configurations that lead to popula¬ 
tions only. In the upper panel of Fig. 1, those contribu¬ 
tions are schematically displayed and grouped according 
to the probability that they occur. 

In contrast to QPT, which is based on 2D-PES, there 
are here fourteen possibilities for the final state instead of 
ten. Thus, twenty independent experiments are needed 
instead of sixteen. This comes at the expense of the dif¬ 
ferent role that the fourth pulse has in each technique, 
namely, heterodyne detection in 2D-PES and the gener¬ 
ation of populations in 2D-FS. However, the signals that 
lead to population of the doubly-excited exciton manifold 
I / ) ( / I from the coherence | / ) ( e | must be summed up 
to the signal that lead to the population | e ) ( e |. In the 
summation, the process |/)(e |—>■ |/)(/| is weighted 
with a factor T, while the process |/)(e| —|e)(e| 
has weight -1. The same procedure applies to the pro¬ 
cess I /)(e'| ^ I / )(/ I and I / )(e'| ^ |e')(e'|. This 
procedure leads to the sixteen independent signals that 
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are needed to reconstruct the sixteen elements of the pro¬ 
cess tensor Xnmu^iiT). 

By using the same toolbox as in the preparation state 
and following closely the notation in Ref. [ 8 ], the total 
signal is given by 




Jei,e2,e3,e4 

E 

p,g,r,sG{e,e'} 


^UJi ^a ;2 ^^3 ^0J4 ^Gi , 62 , 63,04 


(14) 


with 

pp,q,e,e (t T t) = 

{l^pg ■ ei) {Hqg ■ 62 ) ^gp{T) 

X {{fleg ■ es) {fleg ■ e4)'^eg{t) (15) 

^ iXqqqpi'^) ~ ^pq ~ Xeeqp(T)] — (1 — F) 

^ it^fe' ' 63) (R/e' ■ ^4) '^eg {t)Xe' e'qp {T)} 


and 


pp,q,e,e / rp 

^ 61 ,62,63,64 ^ 1 


t) = 


{flpg ■ ei) {/j'qg ■ 62 ) ^gp{T) 

X {{fleg ' G3) {fJ-e'g ' '^e'git)Xe'eqpiT'j 

+ (1 - r) (R/e' ■ 63 ) (/i/e • 64 ) %(t)} 


(16) 


Analogous expressions hold for Pil%lje 3 ,e 4 
-Pei'e 2 es 64 after interchanging e O e'. The 2D-FS 
signalsl 5 'Fs]e(:e 2 %“ 3 %?(T'r,l) in Eq. (14) with (15) and 
(16) are the main result of this article. They allow for 
the reconstruction of the dynamics of excitonic systems 
based on 2D-FS that is an attractive approach to reach 
QPT at the level of single molecules. Remarkably, the 
appealing form of Eqs. (14), (15) and (16) allows for 
an immediate connection with the protocol derived in 
Refs. [ 8 , 9]. Specifically, up to a global minus sign that 
is consistent with previous investigations [16, 20 ], results 
in Refs. [ 8 , 9] are obtained by simply setting F = 2 in 
Eqs. (14), (15) and (16). 

Because the probed sample is an ensemble of isotropi¬ 
cally distributed molecules in solution, an isotropic aver¬ 
age of (/ta ■ ei) (/Xfc • 62 ) (/tc • 63 ) (/td • 64 ) is needed. In 
doing so, standard procedures are followed (see, e.g., 
Chap. 11 in Ref. [27] or Sec. 3.3 in Ref. 13). Specifically, 
the isotropic average is given by 


{{f^a ' Gl) {fJ'b ' ^ 2 ) (/^c * ^ 3 ) {f^d ' G4))iso 


= E 


I ( 4 ) 

61,62,63,64,mi,7712,m3,7714 


7711,m2,m3,m4 


X (/ta • mi) {fib ■ m 2 ) {fic ■ m 3 ) {fid ■ m 4 ) . 


(17) 


where and m.^ denote the polarization of the pulses 
in the laboratory and molecule-fixed frames, respectively. 
Cj = {ea:i,ey 4 ,ea;i} and rui = {mxi, Wyi, ma,,} are the 
components of the polarization vectors and m^, re¬ 
spectively. The isotropically invariant tensor 1*^^^ is given 


by 


I ( 4 ) 

61,62,63,64,mi ,m2,m3,m4 


1 


(<^ 6162^6364 ^6163 



<5eie4 ^6263 ) 


^mim2 ^17137714 
^ 77177714^77127713 


(18) 


The explicit expression for the relevant case of interest 
in the collinear configuration used in 2D-ES, ei = 62 = 
63 = 64 = z, can be found in Appendix C. Thus, after 
isotropically averaging, 

([5Fs]^(:r2%l%7(TT,t))i3o 

— \ ' CP cq c"^ / pp,q,r,s \. (19) 

~ ^ '-^(.di'-^i.c;2'-^l.iJ3^1.i;4\'‘ ei, 62 , 03 ,e 4 /iso- 

p,q,r,sG{e,e'} 


Because in 2D-FS the laser pulses are collinear, it is pos¬ 
sible to set Bj = z at this point. Additionally, it is as¬ 
sumed below that r = 0 and / = 0 so that only the signals 
-Pejfea’.ea ,64 (0,7,0) are Considered. However, following a 
similar procedure as in Ref. [ 8 ], this restriction can be 
relaxed and arbitrary r and t can be considered. Eor the 
sake of generality, the polarization vectors are denoted 
independently by ej and r and t are not set to zero in 
the above expressions. 

The extraction procedure of the matrix elements of x 
from ([,5Fs]“(,’e“,%“%?(r,r,/))iso follows from Eqs. (19), 
(15) and (16). Note that in doing so the sixteen 2D- 
FS signals ([ 5 'Fs]e(;e 2 '',e^% 4 ^('^,7,t))iso, and the sixteen 
auxiliary signals {Pil%l%^e4{'r,T,t))7so can be grouped 
into the sixteen-dimensional vectors ([Sfs](a 7,t))iso 
and (P(r,T, t))iso, respectively. This allows writing 
Eq. (19) as ([SFs]('r,7,t))iso = C(P(r,T,t))iso, where 
the matrix elements of C contains the probabilities 
. Then, the first step in the extraction 
procedure is to invert the matrix C so that the sig¬ 
nals {Pife 2%3 64)130 can be extracted from the measured 
signals ([ 5 'Fs]e(,’S',e 1 "e 4 '‘(T 7, t))iso, i-e., (P(r, T,/))iso = 
C~^([Sfs]('f, 7, t))iso- The second step comprises the ex¬ 
traction of the sixteen elements of the process tensor 
x(T) from the isotropically-averaged version of Eqs. (15) 
and (16). This process can be accomplished by conve¬ 
niently defining a sixteen-dimensional vector x(7), such 
that x(7) = M~^(P(0,T, 0))iso- See the Appendices for 
further details. 


IV. NUMERICAL EXAMPLE 

As a concrete example, consider parameters of rele¬ 
vance in the context of light-harvesting systems [28, 29] . 
Specifically, to compare with previous results [ 8 ] , consider 
1371 = 12881 cm“^, 1372 = 12719 cm“^, J = 120 cm“^, 
^ 2/^1 = 2 and (p = 0.3. The two-waveform toolbox is as¬ 
sumed to have frequencies uj^. = 13480 cm“^ and uj- = 
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12130 cm“^ so that Wi = {w+,a;_}, Vi and pulse width 
cr = 40 fs. To simulate the signals, the spectral density 
Jn{i^) of the local vibrational environments are assumed 
identical and given by Jn{i^) = (Ar.e./wc)a; exp(—w/wc) 
where the cutoff frequency is set as = 120 cm“^ while 
the reorganization energy is chosen as Ar.e. = 30 cm“^. 
These set of parameters are relevant for light-harvesting 
systems and were used in Ref. [8]. 

In the simulations below, an inhomogeneously broad¬ 
ened ensemble of 10^ dimers with diagonal disorder is 
considered. Specihcally, it is assumed that the site ener¬ 
gies w'l and w '2 in the ensemble follow a Gaussian distri¬ 
bution centered at W\ and -072 with standard deviation 
CTinh = 40 cm“^. The dynamics are solved at the level 
of the secular Redfield master equation at room temper¬ 
ature and for T > 3 ct. 

Figure 2 depicts the nonvanishing real parts of 
(^r.z,z,z)iso for a variety of values of the quantum yields 
parameter F. From the functional dependence on F of 
the signals in Eqs. (15) and (16), three cases 

are of interest: (i) For F = 0, the contribution from the 
excited state absorption (ESA) pathways has the same 
sign as the stimulated emission (SE) and ground-state 
bleach (GSB) contributions (see the double-sided Feyn¬ 
man’s diagrams in Fig. 1 or the discussion in Ref. [16]). 
Thus, the amplitude of the signals is the largest possible, 
(ii) For F = 1, the ESA pathways do not contribute to 
the signal and the amplitude of the signals is expected to 
be smaller than in the case of F = 0. (hi) For F = 2, the 
contribution from the ESA pathways has opposite sign 
to the SE and GSB contributions, so that the amplitude 
is expected to be smaller than in the previous case F = 1 
. These expectations are confirmed by simulations in 
Fig. 2. For completeness, the intermediate cases F = 0.5 
and F = 1.5 were also depicted in Fig. 2. 

Based on the signals (T’J'z.z.’Diso obtained above, the 
QPT tensor is reconstructed in Fig. 3. The reconstruc¬ 
tion appears to be insensitive to the value of the quan¬ 
tum yield parameter F. This unexpected result can be 
understood after noticing that the QPT tensor Xnmj//i(T) 
is a characteristic of the singly exited exciton manifold 
and P is a function of the doubly-excited exciton man¬ 
ifold. Thus, QPT of the singly-exited exciton manifold 
by 2D-FS is robust against nonradiative processes of the 
doubly-excited exciton manifold and benefits from the 
quality of the signals discussed above. 

Results in Fig. 3 agree with the secular Redfield tensor 
used to simulate the 2D-FS signals. If the signal were ob¬ 
tained from experimental data, a careful analysis of the 
propagation of errors would be in order. In particular, it 
is necessary to include fluctuations in the laser intensity 
at each time T at which the signals are collected, and 
to pay attention to the stability conditions imposed by 
the invertibility of the matrices C and M [8] . In this arti¬ 
cle, interest was in providing a proof-of-principle for the 
scheme derived above, so that Fig. 3 is aimed to depict 
the type of information that can be extracted from the 
protocol. 



FIG. 2. Nonvanishing T, 0))iso signals for a variety 

of values of F and for Sa < T < 700 fs. Curves with F = 2 
coincides with those extracted in Ref. [8]. 


Specifically, (i) if for a particular photochemical sys¬ 
tem, non-negligible, non-secular terms emerge during the 
reconstruction of the process tensor x from experimen¬ 
tal data, that would imply, e.g, that coherent control 
schemes assisted by the environment [23, 30] may be ap¬ 
plied in that particular system, (ii) If the decay of the 
tensor elements associated to the coherences of the den¬ 
sity matrix, Xnmnm with n ^ m, are non-exponential, 
it may indicate the presence of non-Markovian dynamics 
[31]. The deviation from exponencial decay behavior may 
even be considered as a measure of the non-Markovian 
character of the dynamics- a relevant topic in the con¬ 
text of open quantum systems, (iii) Although in mul- 
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FIG. 3. Nonvanishing elements of the QPT tensor Xnmvij.iT) 
obtained form the depicted in Fig. 2 for a variety 

of values of F and for Sa < T < 700 fs. Curves with F = 2 
coincides with those extracted in Ref. [8]. 

tilevel systems the decay rate of the elements Xnmnm 
cannot be directly associated to the decay rate of the 
coherences (n|/3|m) of the density matrix, the decay rate 
of Xnmnm pi'ovides information about the lifetime of par¬ 
ticular transfer and coherent mechanisms. 


V. DISCUSSION 

Having experimental access to the process tomography 
tensor Xnmnn (t) is fundamental to revealing energy path¬ 
ways in exciton dynamics, and in designing control strate¬ 
gies to increase transport efficiency. Specifically, applica¬ 


tions of QPT to photosynthetic light-harvesting systems 
can (i) rule out certain transfer mechanisms proposed in 
the literature, and (ii) address the question about the 
quantum/classical nature of the energy transport in cer¬ 
tain biological systems from an experimental viewpoint. 
In addressing these issues, a complete analysis of the clas¬ 
sical/quantum correlations encoded in the process tomog¬ 
raphy tensor, as well as an analysis of the main contribut¬ 
ing elements to energy transport is required and will be 
discussed elsewhere. 

A variety of applied and foundational problems can 
be addressed once the process tomography tensor is re¬ 
constructed. From a foundational viewpoint, if the pro¬ 
cess tomography tensor is translated into the phase-space 
representation of quantum mechanics, it reduces to the 
propagator of the Wigner function [32-34]. Based on this 
object, it is possible to experimentally reconstruct signa¬ 
tures of quantum chaos such as scars with sub-Planckian 
resolution [32]. Phase-space resolution below h can be 
achieved here because the process tomography tensor, or 
equivalently the propagator of the Wigner function, is 
not a physical state and therefore, it is not restricted by 
the uncertainty principle [32]. 

In the same way that 2D-PES was extended to study 
chemical exchange to obtain reaction rates under well 
controlled conditions (see, e.g., Chap. 10 in Ref. 13), a 
straightforward extension of QPT is the accurate mea¬ 
surement of concentration of different species in chem¬ 
ical reactions. This has been considered very recently 
in the literature [35]. In this context, interest is in the 
population dynamics Xnnnv(T) of the different chemical 
species, which under Markovian dynamics are in accor¬ 
dance to detailed balance and Onsager’s regression hy¬ 
pothesis [25, 36, 37]. In this respect, because ultrafast 
spectroscopy allows for the study of chemical exchange 
with no need of pressure, temperature, pH nor concen¬ 
tration jumps, it is expected that the proposed approach 
provides experimental evidence for the failure of the On¬ 
sager’s regression hypothesis induced by non-Markovian 
dynamics at the quantum level [25, 36, 37]. 

The QPT scheme introduced here can be readily im¬ 
plemented at the experimental level, and constitutes a 
first step toward the formulation of QPT at the single¬ 
molecule level. Such a scheme would certainly incorpo¬ 
rate quantum aspects of the electromagnetic radiation 
such as the use of time-energy-entangled photons [38]. 
This is already under development in our laboratories. 
Finally, based on present non-linear optical activity spec¬ 
troscopy (see, e.g., Chap. 16 in Ref. 13), by introduc¬ 
ing time-polarization-entangled photons, instead of time- 
energy-entangled photons [38] , single-molecule QPT may 
be extended to study optically active materials at the 
single-molecule level. These materials exhibit unique op¬ 
tical properties, and are constantly finding applications 
in science and industry. 
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Appendix A: Initial State Preparation 


Because the effective initial state is prepared by the first two pulses, it is of second order in A. Thus, 

after applying second order perturbation theory to the time evolution of the system density matrix (see, e.g.. Chap. 5 
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in Ref. [12]), the effective initial state reads 


rei,ei 


{t + T) 


rt2+T 


rt" 


dt" 




(Al) 


where it was assumed that p{—oo) = \g){g\- Symbols in calligraphic font denote superoperators; in particular, 
y = where = \Vi,-]. Assuming that the rotating wave approximation (RWA) holds, and that the 

rephasing signal is synchronously-phase detected at ^ = —4>i + (f >2 + (ps — (pi, the interaction with the electromagnetic 
radiation is conveniently described by 

Cl = -A/i< • eiE{t - (A2) 

C = -AA> • e2E{t - (A3) 

C 3 = -AA> • e^Eit - (A4) 

V 4 = -A/i< • e 4 E{t - , (A5) 


where /i^ = Mp<?b)('7l promotes emissions from the ket and absorptions on the bra, and induces 

the opposite processes. 

In the following, it is considered that the first and the second pulse are well separated, i.e., t = t 2 — ti > Scr. This 
allows for the substitutions y{t") = yit") and y(t") = 


t 2 +T 


. —00 
t" 


P:i::pit + T) = X^Y.\ J 


X J -t')\g){g\(^figg-e2\q){g\E{t" -t2)e 


t2+T 


(A6) 


J dt"x(T) %git2-t'') j dt'%p{t” -t')\g){g\(^fj.pg-ei\q){g\E{t' 


X (^tJ,gg-e4\q){g\E{t''-t2)e 


If the duration of the pulse a is much sorter than the dynamics induced by the environment characterized by Tnm, 
i.e., if (T <C then decoherering contributions can be neglected so that Qgp(ti — t') k. exp [iWpgpti — t')], Gqp{t2 — 
t")Qgp{t'' - 12) ^ exp [-iwqp{t2 - t”)] exp [iwpg{t2 - t")] = exp [-iwqg{t2 - t")], Ggp{t" - < 2 ) ~ exp [-izupg{t" - ^ 2 )] 
and Gggit" — 12) ~ exp [—\Wgg(t" — ^ 2 )]- Moreover, ’^^gpit" — t') ~ ^gpit" — hWgpiti — t')^gp (^2 ~ ^i)- After some 
manipulations, 

PTPep{t + T) ^ -xiT)'^{CP^C;i^^{fj,pg ■ ei){fiqg ■ e2)%p{T) (|q)(p| - 5pq |g)(5l)} (A7) 

pq 


with 


CP. = —— / ds exp[i(wj — Wpg)s]E{s) = — — \/ 2 'Ka'^ exp[—(tjjpg — Wj)^] (A8) 

' 1 J-00 1 

As mentioned in the main text, this effective initial state coincides with the one prepared by 2D-PES in Ref. [8]. 


Appendix B: Final State Detection 

To derive the explicit form of the density operator after the action of the four pulses, it is assumed that the third 
and the fourth pulses are well separated as well. To account for the action of the third pulse and a subsequent period 
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of free evolution, perturbation theory is applied once more, so that the density operator of the system reads 

(t + T + t) = ^ { [^^3 • e 3 )i^eg(t)(ff {t + T)\ g) 

pq 

+ T)| e) - <(/Xe's • e3)i#e'g(t)(e + T)\ e')] \e){g\ 

+ • es)%>g{t){g + T)\g) 

|p--,7(t + T)| e') - CZ,{piea • + T)\e)] \e’){g\ 

+ C'l^iaCM/e • e 3 )i#/e(t)(e |pei\’ei^(^ + T")! e) + C'®3(/X/e/ • e 3 )^#/e W (e' |Pei\’ef (^ + "r) | e) |/)(e| 

+ [<(M/e • B,Wfe'{t){e\&{t + T’)! e') + C:3(M/e' • B,Wfe'{t){e' \ptlTAt + r)| e')] |/)(e'|} , 


(Bl) 


where 


(* \p%l:^: {t + T)\ 3 ) = -Y. CP3 (/Xpg . ei)(M,, 


■ 62 ) 


^9P\ 


('^) {XijqpO^) ^pq^ij^ig) ■ 


pq 


Finally, the fourth pulse prepares the system in the state, 

-a;i,a;2,a;3,a;4 (^ \ \ i\ 

rei,e2,e3,e4 ' ^) 

= ■ e,){f \p:1’ZZ^{t + T + t)\e)- ■ e^){g + T + t)\e)) |e)(e| 

+ {CZfie'f ■ e,){f \pZZZi^ + T + t)\e')- ■ e,){e' + T + t)\g)) \e'){e'\ 

- {ciSpiee'-^Z{f\p%lZZZ + T + t)\e)+CZMe,f-B,){f\ZZZ^^ |/)(/|. 


(B2) 


(B3) 


Each contribution can be easily associated to the double-sided Feynman diagrams in Fig. 1 in the main text. 

Once the state of the system is obtained, the spectroscopy signals [ 5 'Fs]ei\’e 2 ^,<^^e 4 ^('’’j synchronously detected 
at ^ = —(j)i + (j )2 + (j>3 — ^ 4 , follow from the calculation of {A{t + T + t)) = tvAp^Z^ZZiZ + T + t) with A = 
Ei/={e.e',/}B>v \v){v\. Specifically, 


(i(r + r + t)) 


- c:,(Mge • 64) (e I pZZZ{ r + T + t)\g) 

- • 64) (e' \pZZ:ZZ + T + t)\g) 

+ <(/Xe/ • 64)(1 - r)(/ \ pZZ:ZZ + T +1)| e) 

+ C^J/x/e' • 64)(1 - r)(/ \pZZZZ^ + T + t)\e'). 


(B4) 


After replacing the explicit functional form of the density matrix elements {v \ ptlZZZ + T’ + t)| in Eq. (B4) and 
after conveniently collecting terms, Eq. (B4) leads the 2D-FS signals in Eq. (14) that are the main result of this 
article. 


Appendix C: Isotropic Averages 

Before proceeding to the calculation of the isotropic average, it is necessary to express the dipole transition op¬ 
erators in the molecular frame. In doing so, take as reference the transition dipole operator pb^g = PegT^z- Hence, 
pie.'g — Pe'g e'g')i^z “b Me'p Sin(^g/p)in 3 ;, Pfe — A^/e 60 s(^Je)in 2 ; -f //yg sin(^jg)in 3 ; and Pfe' — Pfe' COs(^Jg/)lll 2 ; -f 

Pfe' sin( 0 /e')m 2 ,. The angle between the different transition dipole moments is given by 

tan(0.^) = (Cl) 

Peg • Pvg. 

The isotropically averaged signals can then be written in the compact form 


pPi = MP'^X^P P,q€{e,e'} 


(C2) 
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with 


pee(j.) 

pe'e'(j.) 

(T) 


T, 0))i.o, (0, T, 0))iso, (0, T, 0))i.o, (0, T, 0))i.o , 

{Pi'.ino, T, 0))i,o, (0, T, 0))i,o, T, 0))i,o, (P|,(0, T, 0)) 

(^l'z:=zT(0, T, 0))iso, {Pihr'io. T, 0))iso, {Pi'i^YiO, T, 0))iso, {PlfAViO, T, 0))iso. 

(i"z%^,'zT(0, T, 0))iso, {Pi:ii%^\o, T, 0))iso, (P|;|,'.:=?(0, t, o))i.o, {Pi:i^:/{o, t, o)hso 


(C3) 


and 


X'="(P) = 

X^'^'{T) = 
X^’^'{T) = 


^e,e,e,e(j,)^^e ,e ,e,e (y) ^ ,e.e (j.) ^ ^^e,e 

^^e.e,e,e' CJ^e'.e',e,e' c^^e.e',e,e' (T) ^ ^ 


(C4) 


'll 


'13 


'23 


A 4 

15^^®’ 

'\\l = = Mr2 = ^Ti = 


^ .4 _ n _ T\ ^ 
15 

/lee 

'43 


?2 = - (1 - [C0s(26»/e') + 2] 

= 0 


'32 


'44 


= -TTML {(1 - r) [3 008(0/6) cos( 0 /e/) + sin( 0 /e)sin( 0 /e/)] M/eM/e' + 3 COs(0e's)MesMe'g} 


33 15^®® 


(C5) 


'24 


= -iM 
1 


23’ 




'41 


= “Y^Meg {[Cos(206'g) + 2] g + (1 - T) [cOs(20/6) + 2] , 


'42 


~ [C0s(2^e/^) + 2] 


Ml{^ 

_ "22 

“ g^eg^e' 

Mi/ 

= Mi/ = K 

^i/ 

= = - 

Mi/ 

= --Mi/, 

Mi/ 

= “Y^^e'g 

Mi/ 

2 4 

“ 5^e'g 


le e _ IWie e 

'21 — ''''22 


IVI 31 


^?2^ 


/le e 

'32 


15' 


f\e e _ :|W|ee 

I 34 — IIVI 23 , 


= {[Cos(20e'g) + 2] filg + (1 - T) [cos(20/e' - 9e'g) + 2] 

= Mi/ = Mi/ = 0 

^/e') “1“ COs(^je “1“ ^fe' ^^e'g^] l^fef^fe' H" ^ g)f^egf^e'g} 

de'g) + 2 ] /t/g} 


(C6) 
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/ 2 /I 

Mil = -5 C0s(6>e/g)^^j^//e'g, = - —/ieg/ie'g {3 C0s(6»e/g)/ieg + (1 - T) [cOs( 26 »/e- - 6>e'g) + 2 C0s(6»e/g)]} 

= Mil = MU 

Mli = Uli = = Uli = M“' = M^t' = 0 

M24 = Afgg = Y^MegMe^g {(1 r)[2cOs(^Je ^e^g) “t“ 2 COs(^Je) COs(^^e^ ^e^g)] M/eM/e^} 

M-' = -iM^'/ 


MU = 

= iM^^, 

M 

ee 

31 

= Mr2 

M 

S - M 

ee 

35 


M 

ee 

36 — 

:Mr 

= 0 , 

Mr = 

= iMr, 

Mfi' = 


gA^e'g 

{3 

COs(6>e<, 

'^)^^e'g 

+ (1- 

r) 

[cos(20/e 

-0e' 

g) + 2 COS(0e'g)] 

M^e} 

Mt^' = 

2 

= - - cos 
5 

{Se'g)l^e; 

gf^e'g 











Mfs' = 

= Mli = 

Mf; 


Mfs' = 

0, 

Mtl' 

= i 

M 

ee 

417 


Mt; 

- iMtr 



M-' = 

= M?f, 

M 

ee' 

52 

= Mli 

1 

M^^' = 

m; 

ee' 

54 


Mr; 

M 

rs = 0 

Mr; 

= -iMfi, Mr' 

M-' = 

= Mli = 


== 

Mli = 


= Mr; 

= 

b' 

\ee' 

'68 

= 0 

Mr; = M 

e' e' 

24 ’ 

Mr' = -iM^t', 

= 

= M?!' = 


= 

= 

Myg 

= ivir' 

= 

0, 


IVI78 — 

'^'37 ’ 



M-' = 

= Mii, 

M 

ee' 

82 

= Mf; 

5 










Mii -- 

= Mii = 

Mg? 

— 

M-' = 

0 

'v'85 - 


M 

ee 

45’ 


'''■86 

— IVIgg 




In defining the vector x(T), besides its properties in Eq. (2)-(4), the following relations were instrumental, Xggee{T) — 

1 = -Xeeee{T) - Xe'e'eeiT), Xgge'e'{T) - 1 = -Xeee'e’{T) - Xe'e'e'e'{T) and Xggee'{T) = -Xeeee'{T) - Xe'e'ee'{T). 



